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Abstract 

Motivated by non-linear, non-Gaussian, distributed multi-sensor/agent navigation and tracking appli- 
cations, we propose a multi-rate consensus/fusion based framework for distributed implementation of the 
particle filter (CF/DPF). The CF/DPF framework is based on running localized particle filters to estimate 
the overall state vector at each observation node. Separate fusion filters are designed to consistently 
assimilate the local filtering distributions into the global posterior by compensating for the common past 
information between neighbouring nodes. The CF/DPF offers two distinct advantages over its counter- 
parts. First, the CF/DPF framework is suitable for scenarios where network connectivity is intermittent 
and consensus can not be reached between two consecutive observations. Second, the CF/DPF is not 
limited to the Gaussian approximation for the global posterior density. A third contribution of the paper 
is the derivation of the exact expression for computing the posterior Cramer-Rao lower bound (PCRLB) 
for the distributed architecture based on a recursive procedure involving the local Fisher information 
matrices (FIM) of the distributed estimators. The performance of the CF/DPF algorithm closely follows 
the centralized particle filter approaching the PCRLB at the signal to noise ratios that we tested. 

lEEEkey words: Consensus algorithms. Data fusion. Distributed estimation. Multi-sensor tracking. Non- 
linear systems, and Particle filters. 

I. Introduction 

The paper focuses on distributed estimation and tracking algorithms for non-linear, non-Gaussian, 
data fusion problems in networked systems. Distributed state estimation has been the center of attention 
recently both for Unear |[5l- |[T0l and non-linear systems ifTTI - llBll with widespread applications such 
as autonomous navigation of unmanned aerial vehicles (UAV) |[T3l . localization in robotics |[T5l . track- 
ing/localization in underwater sensor networks |16|, distributed state estimation for power distribution 
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networks |[T4l . and bearings-only target tracking lITTl . A major problem in distributed estimation networks 
is unreliable communication (especially in large and multi-hop networks), which results in communication 
delays and information loss. Referred to as intermittent network connectivity |[33l . |[34l . this issue has been 
investigated broadly in the context of the Kalman filter ll33l . |[34i . Such methods are, however, limited 
to linear systems and have not yet been extended to non-linear systems. The paper addresses this gap. 
Distributed Estimation: Traditionally, state estimation algorithms have been largely centralized with 
participating nodes communicating their raw observations (either directly or indirectly via a multi-hop 
relay) to the fusion centre: a central processing unit responsible for computing the overall estimate. 
Although optimal, such centraUzed approaches are unscalable and susceptible to failure in case the 
fusion centre breaks down. The alternative is to apply distributed estimation algorithms, where: (i) There 
is no fusion center; (ii) The sensor nodes do not require global knowledge of the network topology, and; 
(iii) Each node exchanges data only within its immediate neighbourhood limited to a few local nodes. 
The distributed estimation approaches fall under two main categories: Message passing schemes (TSl, 
|[T9l . where information flows in a sequential, pre-defined manner from a node to one of its neighboring 
nodes via a cyclic path till the entire network is traversed, and; Diffusive schemes lfT4l - |[T7l . |[20l - |[29l . 
where each node communicates its local information by interacting with its immediate neighbors. In 
dynamical environments with intermittent connectivity, where frequent changes in the underlying network 
topology due to mobility, node failure, and link failure are a common practice, diffusive approaches 
significantly improve the robustness at the cost of an increased communications overhead. In diffusive 
schemes, the type of information communicated across the network varies from local observations, 
their UkeUhoods, or some other function of local observations HH, Ell, ||24l, HH, ||29l, ED, to state 
posterior/filtering estimates evaluated at individual nodes fTOll . |[20l . |[22l . |[23l . ll26l . ll27l . Communicating 
state posteriors is advantageous over sharing likelihoods in applications with intermittent connectivity. 
In theory, any loss of information in error prone networks is contained in the following posteriors 
and is, therefore, automatically compensated for as the distributed algorithms iterate. A drawback of 
communicating the local state estimates stems from their correlated nature IB- Channel filters ID and 
their non-linear extensions |[20l (proposed to ensure consistency of the fused estimates by removing 
this redundant past information) associate an additional filter to each communication link and track 
the redundant information between a pair of neighbouring nodes' local estimates. However, channel 
filters are limited to tree-connected topologies and can not be generalized to random/mobile networks. 
Alternatively, consensus-based approaches have been recently introduced to extend distributed estimation 
to arbitrary network topologies with the added advantage that the algorithm is somewhat immune to node 
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and/or communication failures Q, ifTOl . The consensus-basej] distributed Kalman filter ll5l- |[T0l have 
been widely explored for estimation/tracking problems in linear systems with intermittent connectivity 
but there is still a need for developing distributed estimation approaches for non-linear systems. In 
addition, the current non-linear consensus-based distributed approaches suffer from three drawbacks. 
First, the common practice ll2T]| - |[23l |[26l of limiting the global posterior to Gaussian distribution is 
sub-optimal. Second, common past information between neighbouring nodes gets incorporated multiple 
times |[22i . Finally, these approaches ll2ll - |[29l . require the consensus algorithm to converge between 
two successive observations (thus ignoring the intermittent communication connectivity issue in the 
observation framework). The performance of the distributed approaches degrade substantially if consensus 
is not reached within two consecutive observations. 

Motivated by distributed navigation and tracking applications within large networked systems, we pro- 
pose a multi-rate framework for distributed implementation of the particle filter. The proposed framework 
is suitable for scenarios where the network connectivity is intermittent and consensus can not be reached 
between two observations. Below, we summarize the key contributions of the paper. 
1. Fusion filter: The paper proposes a consensus/fusion based distributed implementation of the particle 
filter (CF/DPF) for non-linear systems with non-Gaussian excitation. In addition to the localized particle 
filters, referred to as the local filters, the CF/DPF introduces separate consensus-based filters, referred to 
as the fusion filters (one per sensor node), to derive the global posterior distribution by consistently fusing 
local filtering densities in a distributed fashion. The localized implementation of the particle filter and the 
fusion filter used to achieve consensus are run in parallel, possibly at different rates. Achieving consensus 
between two successive iterations of the local filters is, therefore, no longer a requirement. The CF/DPF 
compensates for the common past information between local estimates based on an optimal non-linear 
Bayesian fusion rule ||35l . The fusion concept used in the CF/DPF is similar to HI and |[20l . where separate 
channel filters (one for each communication link) are deployed to consistently fuse local estimates. Fig. [T] 
compares the proposed CF/DPF framework with channel filter framework and the centralized Architecture. 
In the channel filter framework (Fig. [Jc)), the number of channel filters implemented at each node equals 
the number of connections it has with its neighbouring nodes and, therefore, varies from one node to 
another. These filters are in addition to the localized filters run at nodes. In the CF/DPF (Fig. [IJa)) 
each node only implements one additional filter irrespective of the neighbouring connections. Finally, the 

'Consensus in distributed filtering is the process of establishing a consistent value for some statistics of the state vector across 
the network by interchanging relevant information between the connected neighboring nodes. 
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Fig. 1. (a) The proposed CF/DPF implementation where sensor nodes connect through their fusion filters (one fusion filter per 
node), (b) Centralized implementation where all nodes communicate their local estimates to the fusion center, (c) Distributed 
implementation using channel filters where a separate filter is required for each communication link. In terms of the number of 
extra filters, the CF/DPF falls between the Centralized and channel filters. 

tree-connect network shown in Fig. [TJc) can not be extended to any arbitrary network, for example the 
one shown in Fig. [Ja). The CF/DPF is appUcable to any network configuration. 

2. Modified Fusion filters: In the CF/DPF, the fusion filters can run at a rate different form that of 
the local filters. We further investigate this multi-rate nature of the proposed framework, recognize three 
different scenarios, and describe how the CF/DPF handles each of them. For the worse-case scenario 
with the fusion filters lagging the local filters exponentially, we derive a modified-fusion filter algorithm 
that limits the lag to an affordable delay. 

3. Computing Posterior Cramer-Rao Lower Bound: In order to evaluate the performance of the 
proposed distributed, non-linear framework, we derive the posterior Cramer-Rao lower bound (PCRLB), 
(also referred in literature as the Bayesian CRLB) for the distributed architecture. The current PCRLB 
approaches 12, 13611 . ||37l assume a centralized architecture or a hierarchical architecture Q. The exact 
expression for computing the PCRLB for the distributed architecture is not yet available and only an ap- 
proximate expression @ has recently been derived. The paper derives the exact expression for computing 
the PCRLB for the distributed architecture. Following Tichavsky et al. [2|, we provide a Riccati-type 
recursion that sequentially determines the exact FIM from localized FIMs of the distributed estimator. 

The rest of paper is organized as follows. Section |ll] introduces notation and reviews the centralized 
particle filter as well as the average consensus approaches. The proposed CF/DPF algorithm and the 
fusion filter are described in Section [Till The modified fusion filter is presented in Section |IVl Section IVl 
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derives an expression for computing the PCRLB for a distributed architecture. Section |Vl] illustrates 
the effectiveness of the proposed framework in tracking applications through Monte Carlo simulations. 
Finally, in Section IVIIl we conclude the paper. 

II. BACKGROUND 

We consider a sensor network comprising of N nodes observing a set of Ux state variables x = 
[Xi, . . . , For (1 < / < A^), node / makes a measurement z^^\k) at discrete time instants 
k, (1 < k). The global observation vector is given by z = [z^i)"^, . . . where T denotes 

transposition. The overall state-space representation of the dynamical system is given by 



State Model: 



Observation Model: 
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where ^(•) and C(-) are, respectively, the global uncertainties in the process and observation models. 
Unlike the Kalman filter, the state and observation functions /(•) and g{-) can possibly be nonUnear, 
and vectors ^(■) and C(-) ai^e not necessarily restricted to white Gaussian noise. 
The optimal Bayesian filtering recursion for iteration k is given by 

J P{x{k-l)\z{l:k-l))f{x{k)\x{k-l))dx{k-l) 
P{z{k)\x{k))P{x{k)\z{l:k-l)) 



P{x{k)\z{l:k-1)) 
and P{x{k)\z{l:k)) 



(3) 
(4) 



P{z{k)\z{l:k-1)) 

The particle filter is based on the principle of sequential importance sampling ll38l . |[39ll . a suboptimal 
technique for implementing recursive Bayesian estimation (Eqs. ([3]) and ^) through Monte Carlo sim- 
ulations. The basic idea behind the particle filters is that the posterior distribution P{x{0:k)\z{l:k)) is 
represented by a collection of weighted random particles {Xj(/i:)}^°]^ derived from a proposal distribution 
q{x{0:k)\z{l:k)) with normalized weights Wi{k) = ^(^^^o'^kllzivk)) ' (1 < * < ^s), associated with 
the vector particles. The particle filter implements the filtering recursions approximately by propagating 
the weighted particles, (1 < i < Ns), using the following recursions at iteration k. 



Time Update: 
Observation Update: 



^i{k) 

Wi{k) oc Wi{k 



q{Xi{k)\X,{0:k-l),z{l:k) 

p(z{k)\Mk))p{Xi{k)\Mk-r 



q{Mk)\XiiO:k-l),zil:k) 
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where oc stands for the proportional sign and the proposal distribution satisfies the following factorization 

'x{0:k)\z{l:k)) = q(x{0:k-l)\z{l:k-l))q(x{k)\x{0:k-l),z{l:k)). (7) 



The accuracy of this importance sampling approximation depends on how close the proposal distribution 
is to the true posterior distribution. The optimal choice B3l for the proposal distribution that minimizes 
the variance of importance weights is the filtering density conditioned upon x{0 : k — 1) and z{k), i.e., 



qi^x{k)\x{0:k-l),z{l:k)j = P(x{k)\x{0:k-1), z{k)y (8) 

Because of the difficulty in sampling Eq. dS), a common choice B3l for the proposal distribution is 
the transition density, P{x{k)\x{k — 1)), referred to as the sampling importance resampling (SIR) filter, 
where the weights are pointwise evaluation of the likelihood function at the particle values, i.e., Wi{k) a 

Wi{k-i)Piz{k)\Mk)). 

B. Average Consensus Algorithms 

The average consensus algorithm [91, ifTOl considered in the manuscript is represented by 

+ 1) = Uuit)xPit) + Ui,it)X^^Ht), (9) 

where Xc\t) is the consensus state variable(s) at node for (1 < / < N), t is the consensus time index 
that is different from the filtering time index k, and represents the set of neighbouring nodes for 
node The convergence properties of the average consensus algorithms are reviewed in |40|. Please refer 
to 191, ITOll for further details on consensus algorithms. 

III. The CF/DPF Implementation 

The CF/DPF implementation runs two localized particle filters at each sensor node as shown in Fig. [T] 
The first filter, referred to as the local filter, comes from the distributed implementation of the particle filter 
described in Section ITlI-AI and is based only on the local observations : k). The CF/DPF introduces 
a second particle filter at each node, referred to as the fusion filter, which estimates the global posterior 
distribution P{x{0:k)\z{l:k)) from local distributions P{x{k)\z^^\l:k)) and P{x{k)\z'^^\l:k-1)). 

A. Distributed Configuration and Local Filters 

Our distributed implementation is based on the following model 

x{k) = f{x{k-i))+m (10) 

z^^\k) = g^^\x{k)) + C^^^\k), (11) 
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for sensor nodes (1 < / < N). In other words, the entire state vector x{k) is estimated by running 
one locaUzed particle filter at each node. These filters, referred to as the local filters, come from the 
distributed implementation of the particle filter and are based only on local observations z^'-^l : k). In 
addition to updating the particles and their associated weights, the local filter at node I provides estimates 
of the local prediction distribution P{x{k)\z^'-\l:k — l)) from the particles as explained below. 

Computation and Sampling of the Prediction Distribution: From the Chapman-Kolmogorov equation 
(Eq. (|3]l), a sample based approximation of the prediction density P{x{k)\z^''\l:k — 1)) is expressed as 

P (a;(fc)|2W(l:/b-l)) =Y^wP{k-l)P (x{k)\xf\k-l)^ , (12) 

i=l 

which is a continuous mixture. To generate random particles from such a mixture density, a new sample 
X.f\k\k — 1) is generated from its corresponding mixture P{x{k)\xf\k — 1)) in Eq. (fT2l) . Its weight 
wP{k-l) is the same as the corresponding weight for 'Kf\k — 1). The prediction density is given by 

P(x{k)\z^^\l:k-1)'^ = Y^wP{k-l)6(^x{k)-Xf\k\k-l)') . 

i=l 

Once the random samples are generated, the minimum mean square error estimates (MMSE) of the 
parameters can be computed. 

B. Fusion Filter 

The CF/DPF introduces a second particle filter at each node, referred to as the fusion filter, which 
computes an estimate of the global posterior distribution P{x{0:k)\z{l:k)). Being a particle filter itself, 
implementation of the fusion filter requires the proposal distribution and the weight update equation. 
Theorem [T] If35l expresses the global posterior distribution in terms of the local filtering densities, which 
is used for updating the weights of the fusion filter. The selection of the proposal distribution is explained 
later in Section III-E. In the following discussion, the fusion filter's particles and their associated weights 
at node I are denoted by {xf'^^\k),wj;^'^\k)}f^l. 

Theorem 1. Assuming that the observations conditioned on the state variables and made at node I are 
independent of those made at node j, (j ^ I), the global posterior distribution for a N— sensor network is 

P[x[^:k)\z{\:k)\ a— ^xP a;(0:A:)|^(l:A:-l) , (13) 

^ ' n^Ii^(^(^)l^«(l:A:-l)) ^ ' 

where the last term can be factorized as follows 

P{x{<d:k)\z{l:k-l)^ = P[x{k)\x{k-l)^p[x{id:k-l)\z{l:k-l)y (14) 
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The proof of Theorem [T] is included in Appendix [A] Note that the optimal distributed protocol defined 
in Eq. (fT3] ) consists of three terms: (i) Product of the local filtering distribution YiiLi P{x{k)\z^^^ {l:k)) 



which is again only based on the local observations and represent the common information between 
neighboring nodes, and; (iii) Global prediction density P {x{0:k)\z{l:k — l)) based on Eq. (fT4l ). The 
fusion rule, therefore, requires consensus algorithms to be run for terms (i) and (ii). The proposed 
CF/DPF computes the two terms separately (as described later) by running two consensus algorithms at 
each iteration of the fusion filter. An alternative is to compute the ratio (i.e., proportional to the local 
likelihood) at each node and run one consensus algorithm for computing the ratio term. In the CF/DPF, 
we propose to estimate the numerator and denominator of Eq. ([T3] ) separately because maintaining the 
local filtering and prediction distributions is advantageous in networks with intermittent connectivity as 
it allows the CF/DPF to recover from loss of information due to delays in convergence. Maintaining the 
likelihood prevents the recovery of the CF/DPF in such cases. 

C. Weight Update Equation 

Assume that the local filters have reached steady state at iteration k, i.e., the local filter's computation is 
completed up to and including time iteration k where a particle filter based estimate of the local filtering 
distribution is available. The weight update equation for the fusion filter is given by 



The CF/DPF is derived based on the global posterior P{x{0 : k)\z{l : k)) which is the standard approach 
in the particle filter literature [38]. Further, we are only interested in a filtered estimate of the state 
variables P{x{k)\z{l : k)) at each iteration. Following |[38l we, therefore, approximate q{x{k)\x{\ : 
k—\),z{l:k)) = q{x{k)\x{k—\) , z{k)) . The proposal density is then dependent only on a;(A;) and z(A;). 
In such a scenario, one can discard the history of the particles X-''^^^(0: A;— 2) at previous iterations [381. 
Substituting Eq. ([T4l ) in Eq. ([T3] ) and using the result together with Eq. ((T]) in Eq. ([TSl) . the weight update 
equation is given by 



which depends on local observations; (ii) Product of local prediction densities Hi^i P{x{k)\z^^^ (1 : k—1)), 




(15) 



w!>^'^^\k) oiwf'^^\k-i) 



Y^^,p{^r\k)\^'Hi:k)) p {i6r\k)\i6r\k-i)) 



(16) 



niIiP(xr'^(A:)|^W(l:A;-l)) g(xp^)(fc)|Xp^)(fc-l),z(fc) 



where 




(17) 
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Observe that only the first fraction in Eq. (fT6l ) requires all nodes to participate. Given the weights 
wj^^'^^\k — l) from the previous iteration, Eq. (fT6l) requires the following distributions 

N N 

Y[p(xf'^\k)\z^^\l:k)') and J]P (xp''^(A;)|zW(l:A;-l)) . (18) 
1=1 1=1 

The numerator of the second fraction in Eq. (fT6l) requires the transitional distribution P{x{k)\x{k — 1)), 
which is known from the state model. Its denominator requires the proposal distribution q{x{k)\x{k — 
l),z{k)). Below, we show how the three terms (Eq. (fTSl ) and the proposal distribution) are determined. 

D. Distributed Computation of Product Densities 

Distributions in Eq. (fTSl ) are not determined by transferring the whole particle vectors and their 
associated weights between the neighboring nodes due to an impractically large number of information 
transfers. Further, the localized posteriors are represented as a Dirac mixture in the particle filter. Two 
separate Dirac mixtures may not have the same support and their multiplication could possibly be zero. If 
not, the product may not represent the true product density accurately. In order to tackle these problems, a 
transformation is required on the Dirac function particle representations by converting them to continuous 
distributions prior to communication and fusion. Gaussian distributions |[T3l . ifTSl . lfT6ll . |[22l . ll23l . 
|[26l . grid-based techniques IHI, Gaussian Mixture Model (GMM) |[T9l and Parzen representations |[20l 
are different parametric continuous distributions used in the context of the distributed particle filter 
implementations. The channel filter framework |[20l fuses only two local distributions, therefore, the 
local pdfs can be modeled [20] with such complex distributions. Incorporating these distributions in the 
CF/DPF framework is, however, not a trivial task because the CF/DPF computes the product of N local 
distributions. The use of a complex distribution like GMM is, therefore, computationally prohibitive. 

In order to tackle this problem, we approximate the product terms in Eq. (fT6l ) with Gaussian distribution 
which results in local filtering and prediction densities to be normally distributed as follows 

p(a;(fc)|zW(l:fc)) oc pW(/fc)) and P (a;(A:)|zW(l:A: - 1)) oc (i/(')(A:), P«(A:)) , (19) 

where ^^^\k) and P^^\k) are, respectively, the mean and covariance of local particles at node / during 
the filtering step of iteration k. Similarly, v^^\k) and R^^\k) are, respectively, the mean and covariance 
of local particles at node / during the prediction step. It should be noted that we only approximate the 
product density for updating the weights with a Gaussian distribution and the global posterior distribution 
is not restricted to be Gaussian. The local statistics at node I are computed as 

fi(^\k)=Y^wP {k)xf\k) and P^^\k)=Y^wPik)(xf\k)-ti'-^\k)^ (xf\k)-ti'-^\k)y (20) 

1=1 i=l 
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Reference PTI shows that the product of N multivariate normal distributions is also normal, i.e., 

N N 

1=1 1=1 

where C is a normalization term (Reference BTI includes the proof). Parameters ^{k) and P{k) are 

^ -1 ^ -1 

P{k) = {^(p^^Hk)) and fx{k) = P{k)x^(p'-^\k)) ii^^\k). (22) 



1 = 1 ^ V ■' 1 = 1 



XW(0) a;«(0) 

Similarly, the product of local prediction densities (Term (fTSl l) is modeled with a Gaussian density 

M{x{k);v{k), R{k)), where the parameters v{k) and R{k) are computed as follows 

^ -1 ^ -1 

R{k) = (^R^^\k)^ and v{k) = R{k) x ^ (^R^^\k)') v^^\k). (23) 

The parameters of the product distributions only involves average quantities and can be provided using 
average consensus algorithms as follows: 

(i) For (1 < / < N), node / initializes its consensus states to X^^{0) = (p(')(fc))"\ x^^{0) = 
{P(^\k)y^n(^\k), X^^{0) = and a;2(0) = (i?W(A;))-it^«(A:), and then Eq. © 
is used to reach consensus with X^^ (t) used instead of Xc\t) in Eq. ^ for the first consensus 
run. Similarly, a;^2 (0 is used instead of xP (t) for the second run and so on. 

(ii) Once consensus is reached, parameters /i*^'^ (A;) and P^^^ (k) are computed as follows 

P{k) = 1/iV X ^lim and /x(A;) = m^n | (x^f (t)) a;g(t)| (24) 

R{k) = yN xlimU xS{t)y'\ and v{k) = limUx^ (t))'' x x<il{t)\ . (25) 



Based on aforementioned approximation, the weight update equation of the fusion filter (Eq. (1161) ) is 

V(xP''(k);«(k),fi(fc))<,{xf'^''>(fc)|xP"(k-l),z(k)) 
Eq. (|26l ) requires the proposal distribution q{x{k)\x{k — l), z{k)), which is discussed next. 

E. Proposal Distribution 

In this section, we describe three different proposal distributions in CF/DPF. 
1. SIR Fusion Filter: The most common strategy is to sample from the probabilistic model of the state 
evolution, i.e., to use transitional density P{x{k)\x{k+1)) as proposal distribution. The simplified weight 
update equation for the SIR fusion filter is obtained from Eq. (l26l ) as follows 



M{xf^''\k);v{k),R{k)) 
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This SIR fusion filter fails if a new measurement appears in the tail of the transitional distribution or 
when the likelihood is too peaked in comparison with the transitional density. 

2. Product Density as Proposal Distribution: We are free to choose any proposal distribution that 
appropriately considers the effect of new observations and is close to the global posterior distribution. 
The product of local filtering densities is a reasonable approximation of the global posterior density as 
such a good candidate for the proposal distribution, i.e., 

N 

qix{k)\x{k-l),z{l:k)) ^YIp (x{k)\z^^'^{l:k)^ , (28) 

1=1 

which means that we generate particles {X^ ' {k)}-_!^ are generated from J\f {^{k) , P {k)) . In such a 
scenario, the weight update equation (Eq. (l26l )) simplifies to 

W^-^i,)^w!-^(k-n'^;^Jf>^'''^'-'^\ (29) 

Next we justify that the product term is a good choice and a near-optimal approximation of the optimal 
proposal distribution (Eq. dSjl). Assume at iteration k, node /, for (1 < / < A^) computes an unbiased 
local estimate x^^\k) of the state variables x{k) from its particle -based representation of the filtering 
distribution with the corresponding error and error covariance denoted by l^^l\k) = x{k) - x^^\k) and 
P'^^\k). When the estimation error /S.^) {k) and IS.x\k), for (1 < i,j < N) and i ^ j are uncorrelated, 
the optimal fusion of N unbiased local estimates x^^\k) in linear minimum variance scene is shown Il42l 
to be given by 

^ -1 ^-1^-1 

P(A;) = (^p(0(A;)) and x{k) = (^P^^\k)^ x J]] (pW(A;)) i«(A:). (30) 

1=1 1=1 1=1 

where x{k) is the overall estimate obtained from P{x{k)\z{l : k)) with error covariance P{k). Eq. ( [30l ) is 
the same as Eq. (l22l ). which describes the statistics of the product of N normally distributed densities. The 
optimal proposal distribution is also a filtering density 113811 . therefore, the proposal distribution defined 
in Eq. ( |28] ) is a good choice that simplifies the update equation of the fusion filter. Further, Eq. ( |28l ) 
is a reasonable approximation of the optimal proposal distribution. From the framework of unscented 
Kalman filter and unscented particle filter, it is well known ll43l that approximating distributions will 
be advantageous over approximating non-linear functions. The drawback with this proposal density is 
the impractical assumption that the local estimates are uncorrelated. We improve the performance of the 
fusion filter using a better approximation of the optimal proposal distribution, which is described next. 

3. Gaussian Approximation of The Optimal Proposal Distribution: We consider the optimal solution 
to the fusion protocol (Eq. (Tdh ) when local filtering densities are normally distributed. In such a case. 
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P{x{0:k)\z{l:k — 1)) is also normally distributed |[35l with mean x^'-'^^\k) and covariance p(''^^)(fc) 

, N N 

p(W-\A;) = +^pOy\k)-Y,R^'^'\k) (31) 



^ ■' '• 



-1 ^ ^ -1 

= p(W-^(yt)[('p(0(A;)') v^^\k) + ^P'^^^'\k)tx^^\k)-^(R^^\k)) v^^\k)].02) 



0=1 i=i 



The four terms x^^^{oo), a; [.2(00), x^^{oo), and a;[.2(oo) are already computed and available at local 



nodes as part of computing the product terms. Fusion rules in Eqs. (|3T| ) and (|32l ) are obtained based 
on the track fusion without feedback |35|. In such a scenario, particles ^(/c) are drawn from 
M{x^^'^^\k),P^^''^\k)) and the weight update equation (Eq. (|29l)) is given by 

V{X<"'(fc);„(fc),fl(t))A-(X<'-">(t);x<lf,P^«) 
The various steps of the fusion filter are outlined in Algorithm [U The filtering step of the CD/DPF is 
based on running the localized filters at each node followed by the fusion filter, which computes the 
global posterior density by running consensus algorithm across the network. At the completion of the 
consensus step, all nodes have the same global posterior available. As a side note to our discussion, we 
note that the CF/DPF does not incorporate any feedback from the fusion filters to the localized filters to 
provide sufficient time for the fusion filter to converge. The main advantage of the feedback is to reduce 
the error of the local filters which will be considered as future work. Finally, a possible future extension 
of the CF/DPF is to use non-parametric models, e.g., support vector machines (SVM) ifTTI . ll32l . instead, 
for approximating the product terms. An important task in CF/DPF is to assure that the localized and 
fusion filters do not lose synchronization. This issue is addressed in Section [TVl 

F. Computational complexity 

In this section, we provide a rough comparison of the computational complexity of the CF/DPF 
versus that of the centralized implementation. Because of the non-linear dynamics of the particle filter, 
it is somewhat difficult to derive a generalized expression for its computational complexity. There 
are steps that can not be easily evaluated in the complexity computation of the particle filter such 
as the cost of evaluating a non-Unear function (as is the case for the state and observation mod- 
els) fl4l . In order to provide a rough comparison, we consider below a simplified linear state model 
with Gaussian excitation and uncorrected Gaussian observations. Following the approach proposed 
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Algorithm 1 Fusion FiLTER({Xf ' ^(A: - 1), Wf '{k - 
Input: {xf'^^''(A; — l),wf'^^\k — - Fusion filter's particles and associated weights. 

Output: {X) ' {k), Wl ' Fusion filter's updated particles and associated weights. 

1: for / = 1 : N, do 

(A:), P(0 (A;), t^W (k) , (k)) = LocalFilter ({Xf ^ - 1) , {k - l)}^^^ , (A:)) 

2: end for 

3: DoFusion({Ai(')(A;),pW(A;)}/^^) computes {fi^^'^^\k), P'^^'^^\k)} for numerator of Eq. 
4: DoFusion({t;(')(A;),fiW(A;)}/IJ computes {v^^^'^^^k), R^^'^^^k)} for denominator of (O. 
5: for i = 1 : N, do 

{(I FF) 1 ^l,¥F 

X • ' (A;) > by sampling proposal distribution defined in Section III-E. 
J i=l 

• Compute weig hts W^^'^^\k) using Eq. (|29]l. 
6: end for 

7: If degeneracy observed ({xf '^^^(A;), H^^''^^^(A;)}£ffJ = Resample({xf'^^^(A;), H^^''^^^(A;)}£f^J. 



in ll44l . the computational complexity of two implementations of the particle filter is expressed in terms 
of flops, where a flop is defined as addition, subtraction, multiplication or division of two floating 
point numbers. The computational complexity of the centralized particle filter for A^-node network 
with Ns particles is of O {{nl + N)Ns). The CF/DPF runs the local filter at each observation node 
which is similar in complexity to the centralized particle filter except that the observation (target's 
bearing at each node) is a scalar. Setting = 1, the computational complexity of the local filter is 
of O (n^A^Lp) per node, where A'lf is the number of particles used by the local filter. There are two 
additional components in the CF/DPF: (i) The fusion filter which has a complexity of □(n^.A^Fp) per 
node where A'^pp is the number of particles used by the fusion filter, and; (ii) The CF/DPF introduces 
an additional consensus step which has a computational complexity of 0{n'^AgNc{U)). The associated 
convergence time Nc{U) = 1/ log(l/rasym(t^)) provides the asymptotic number of consensus iterations 
required for the error to decrease by the factor of 1/e and is expressed in terms of the asymptotic 
convergence rate rasym(C^)- Based on |40|, Nc{U) = — l/max2<i<iv log(|Aj(L'")|), where \i{U) is 
the eigenvalue of the consensus matrix U. The overall computational complexity of the CF/DPF is, 
therefore, given by max |0(A^n^(A'^LP + NFF)),0{n'^AgNc{U))} compared to the computational com- 
plexity O ((n^ + N)Ns) of the centralized implementation. Since the computational complexity of the 
two implementations involve different variables, it is difficult to compare them subjectively. In our 
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simulations (explained in Section IVll ). the value of the variables are: rix = 4:, N = 20, Ng = 10, 000, 
Nlf = NpF = 500, and Nc{U) = 8 for the network in Fig. [Sja) which results in the following rough 
computational counts for the two implementations: Centralized implementation: 3.6 x 10^, and CF/DPF: 
3.4 X 10^. This means that the two implementations have roughly the same computational complexity 
for the BOT simulation of interest to us. We also note that the computational burden is distributed 
evenly across the nodes in the CF/DPF, while the fusion center performs most of the computations in the 
centralized particle filter. This places an additional power energy constraint on the fusion center causing 
the system to fail if the power in the fusion center drains out. 

IV. Modified Fusion Filter 

In the CF/DPF, the local filters and the fusion filters can run out of synchronization due to intermittent 
network connectivity. The local filters are confined to their sensor node and unaffected by loss of 
connectivity. The fusion filters, on the other hand, run consensus algorithms. The convergence of these 
consensus algorithms is delayed if the communication bandwidth reduced. In this section we develop 
ways of dealing with such intermittent connectivity issues. First, let us introduce the notation. We assume 
that the observations arrive at constant time intervals of AT. Each iteration of the local filters is performed 
within this interval, which we will refer to as the local filter's estimation interval. The duration (the fusion 
filters's estimation interval) of the update cycle of the fusion filter is denoted by Tc. Fig. |2] illustrates three 
scenarios dealing with different fusion filter's estimation intervals. Fig. ^a.) is the ideal scenario where 
Tc < AT and the fusion filter's consensus step converges before the new iteration of the local filter. 
In such a scenario, the local and fusion filters stay synchronized. Fig. |2tb) shows the second scenario 
when the convergence rate of the fusion filter varies according to the network connectivity. Under regular 
connectivity Tc < AT and limited connectivity loesses, the fusion filters will manage to catch up with the 
localized filters in due time. Fig. I^c) considers a more problematic scenario when Tc > AT. Even with 
ideal connectivity, the fusion filter will continue to lag the localized filters with no hope of its catching 
up. The bottom two timing diagrams in Fig. IJ^c) refer to this scenario with Tc = 2AT. As illustrated, the 
lag between the fusion filter and the localized filters grows exponentially with time in this scenario. An 
improvement to the fusion filter is suggested in the top timing diagram of Fig. I^c), where the fusion filter 
uses the most recent local filtering density of the localized filters. This allows the fusion filter to catch 
up with the localized filter even for cases Tc > AT. Such a modified fusion implementation requires an 
updated fusion rule for the global posterior density, which is considered next. 

At iteration A; + m, we assume that node /, for (1 < / < N), has a particle-based approximation of 
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Fig. 2. Multi-rate implementation of the local and fusion filters, (a) The ideal scenario where the fusion filter's consensus 
step converges before the new iteration of the local filter, (b) The convergence rate of the fusion filter varies according to the 
network connectivity, (c) The lag between the fusion filter and the local filter grows exponentially. 



the local filtering distributions P{x{k + m)\z^'-\l : k + m)), while its fusion filter has a particle-based 
approximation of the global posterior distribution P{x{0:k)\z{l:k)) for iteration k. In other words, the 
fusion filters are lagging the localized filters by m iterations. In the conventional fusion filter the statistics 
of P(a;(/c+l)jz(')(l:fc+l)), for (1 < / < iV) are used in the next consensus step of the fusion filter which 
then computes the global posterior P{x{0:k+l)\z{l:k+l)) based on Theorem[I] The modified fusion 
filter uses the most recent local filtering distributions P{x{k+m)\z^^\l : k+m)) according to Theorem|2l 

Theorem 2. Conditioned on the state variables, assume that the observations made at node I are 
independent of the observations made at node j, (j ^ I). The global posterior distribution for a N— sensor 
network at iteration k + m is then given by 

P {x{0:k + m)\z{l:k+m)) oc 

TT i^-^ij TT pU(k')\xik'-l))xP{x{0:k)\z{l:k)).{34) 

The proof of Theorem [2] is included in Appendix |B] In the consensus step of the modified fusion 
filter, two average consensus algorithms are used to compute Hi^i rifc-tlfc+i • ^')) ^^'^ 
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N k+m N 

n n P(x{k')\z^^\l:k')^ o^]jM{tx^^\k + l-.k + m),P^^\k+l:k + m)) (35) 

1=1 k'=k+l 1=1 
N k+m. N 

and Yl n P [x{k')\z^^\l:k'-l)^ (xYlj\f{v^^'>{k+l:k+m),R^^'^{k+l:k+m)), (36) 

l=lk'=k+l 1=1 

instead of computing Ylj^-^^ P{x{k)\z'^''\l: k)) and fliLi P{^{k)\z^''\l'- k — 1)) as was the case for the 
conventional fusion filter. The modified fusion filter starts with a set of particles xf^^^'''\k), W^^^^'''\k) 
approximating P{x{0 : k)\z{l : k)) and generates updated particles xf^^^'''\k+m),wj^^^^'''\k+m) for 
P{x{0:k+m)\z{l:k+m)) using the following weight update equation 

^(i,MPP)^^^^^^^a,MPP)^^^ X -—— ^ 1 , (37) 

J\f{Xf' > {k+my,v{k+l: k+m),R{k+l: k+m)) 
which is obtained directly from Eq. (l34l) . Note that the normal approximation in Eqs. (|35])-(|37]) are similar 
to the ones used in the conventional fusion filter. Furthermore, we note that the modification requires 
prediction of the particles from iteration k all the way to fc+m in order to evaluate the second term on 
the right hand side of Eq. (|37] ). Algorithm |2] outlines this step and summarizes the modified fusion filter. 



V. The Posterior Cramer-Rao Lower Bound 

Considering the non-linear filtering problem modeled in Eqs. ([T]) and ^ and the posterior distribution 
(Eq. ([T3] )) used in developing the CF/DPF, the section computes the Posterior Cramer-Rao lower bound 
(PCRLB) for the distributed architecture. We note that the PCRLB is independent of the estimation mech- 
anism and the bound should be the same for both centralized and distributed architectures. The question is 
whether the centralized expressions for computing the PCRLB are applicable to compute the PCRLB for 
other topologies, i.e., the hierarchical and distributed (decentralized) architectures. Reference |3] considers 
a hierarchical architecture with a central fusion center and shows that the centralized expressions can be 
used directly for the hierarchical case. The same authors argue in im that the centralized expressions are 
no longer applicable for distributed/decentralized architectures. The exact expression for computing the 
PCRLB for the distributed architecture is not yet available and only an approximate expression Q has 
recently been derived. In this section, we derive the exact expression for computing the PCRLB for the 
distributed topology. We note that our result is not restricted to the particle filter or the CF/DPF but is 
also applicable to any other distributed estimation approach. 
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Algorithm 2 Modified Fusion Filter 
Input: {xf'^^^\k), wj;'''^^^\k)}^J^'^^^ - Fusion filter's particles and associated weights. 
Output: {xf'^^^\k+m),W^'''^^^\k+m)}^^ updated particles and associated weights. 
I: for k' = k + 1 : k + m, do 

Af{fi^^\k'),P^^\k')) = SaveGaussian({xf)(A;'),w/'^(A;')}^\) 
M{v(^\k'),R^^\k')) = SaveGaussian({xf^(A;' + l|fc'),W^^'^(A;')}£i) 
2: end for 

3: AA(/xW(A: + l:fc+m),P«(/t+l:A;+m)) = SaveGaussian(^nfct™+i AA(/xW(A:'), -P^'H't')))- 
4: A/'(t^«(fc+l:A;+m),i?W(A;+l:A;+m)) = SaveGaussian^ntt^+i -A/'(t^('H't'), -R^'H^'))) ■ 
5: {/x('W(/c+l:/c+m),p('''^PP)(fc+l:fc+m)} = DoFusion({/x(')(/t + l:A:+m),P(')(/c + l:/t + m 
6: {v^^'^^^\k+l:k+m),R^^'^^^\k+l:k+m)} = DoFusion{{v^^\k + l:k+m),R^^\k+h 
7: for i = 1 : A'^FF, do 
8: for fc' = : do 

xf'^PP)(A:') ~ P{x{k')\xf''^^^\k'-l)). 
9: end for 

Xf'^^^\k+m) M{fi^^'^^\k + l:k + m),P^^'^^\k + l:k + m)). 
Compute weig hts W}^'^^^\k+m) using Eq. ([37]). 
10: end for 



The PCRLB inequahty states that the mean square error (MSE) of the estimate x{0:k) of the state 
variables x{0:k) is lower bounded by 

E{{x{0 : k)-x{0 : k)){x{0 : k)-x{0 : k)f}>[J{x{0 : k))]~^ , (38) 

where E is the expectation operator. Matrix J{x{0 : k)) is referred to as the FIM |2| derived from the 
joint probability density function (PDF) P{x{0:k),z{l:k)). Let V and A, respectively, be operators of 
the first and second order partial derivatives given by Va.(fc) = [oj^Xk)' • • • ' dx ^ (k) ] ^^'^ ^xlfc-i) ~ 
^ x{k-i)^]c(k)- form of the Fisher information matrix J{x{0:k)) is defined as Q 

J{xi0:k)) =E{-Al^^:'l]logPixiO:k),zil:k))}. (39) 

An alternative expression for the FIM can be derived from the factorization P{x{0 : k),z{l : k)) = 
P{x{Q:k)\z{l:k)) x P{z{l:k)). Since P{z{l:k)) is independent of the state variables, we have 

J{x{^:k)) =E{ -A^[;;^jlogP(aj(0:A;)l2(l:A:))}. (40) 
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We now describe the centralized sequential formulation of the FIM. 

A. Centralized computation of the PCRLB 

Decomposing x{0:k) as x in J(a;(0:fc)), Eq. (l40l ) simplifies to 



J{x(0:k)) 



A^\k) Ai2(A;) 
A^^k) A22(A;) 



1 X 



. £C(0 



fc-1) 
fc-1) 



x{0 
x(k) 



k—1) ^x(k) 



x\k) 
x{k) 



logP{x{0:k)\zil:k)) \, 



provided that the expectations and derivatives exist. The information submatrix J{x{k)) for estimating 
x{h) is given by the inverse of the {hx x Ux) right-lower block of \j{x{{) : A;))] ^. The information 
submatrix is computed using the matrix inversion lemma lO and given by 

J{x{k)) = A22(A:) - A^^{k){A^^{k)y^A^^{k). (41) 

Proposition [T] 121 derives J{x{k)) recursively without manipulating the large {kux x knx) matrix A^^{k). 
The initial condition is given by J(a;(0)) = E{— Ai^lUj log P(a;(0))}. 

Proposition 1. The sequence {j[x{k)^ } of local posterior information sub-matrices for estimating state 
vectors x{k) at node I, for (1 < I < N), obeys the following recursion 

(42) 

(43) 
(44) 



J{x{k + 1)) = D'^\k) - D'^\k)(^j{x{k)) + D^\k)j D^\k) 
where D^^{k) = E{ - A^j^j log P {x{k + l)\x{k))} 



x(k) 

T TTT, r A X{k + 1) 



D^\k) = [D^\k)Y=W.{- /\ 

c{A;+l) 
Kfc+l) 



x{k) 



\ogP{x{k + l)\x{k))] 



=E{ _ A^[^+;jlogP(^(fc + l)|x(A;))} +E{ - A^J^+;]logP(z^ }, (45) 



J{z{k+1)) 

The proof of Proposition [T] is given in |f2!|. Conditioned on the state variables, the observations made 
at different nodes are independent ,therefore, J{z{k + 1)) |[36l in Eq. (1451 ) is simplifies to 



N 



N 



J{z{k + l)) = J]j(^«(fc + 1)) = J]E{-A^[^+;]logP(z«(fe+l)|^(fc+l))} 



1=1 



1=1 



In other words, the expression for J{x{k + 1)) (Eq. (l42l )) requires distributed information (sensor 
measurement) only for computing J{z{k + 1)). Other terms in Eq. (|42]) can be computed locally. In 
the next section, we derive the distributed PCRLB. 

B. Distributed computation of the PCRLB 

In the sequel, j(') (a;(0 : k)^ , for (1 < Z < A^), denotes the local FIM corresponding to the local estimate 
of x{{):k) derived from the local posterior density P(a;(0 : A;)!^^') (1 : A;)). Similarly, j(') (a;(0 : A+1|A;)) 



September 6, 2012 



DRAFT 



19 



denotes the local FIM corresponding to the local prediction estimate of x{0:k+l) derived from the local 
prediction density P{x{0:k + l),z^'-^l:k)). The expressions for J^''^ [x{0 : k)) and (a;(0 : fc + l|fc)) 
are similar in nature to Eq. (|4T]) except that the posterior density P{x{0 : k)\z{l : k)) is replaced by 
their corresponding local posteriors. The local FIM j(') (a3(A;)) is given by the inverse of the {ux x n^,) 
right-lower block of [j^'^ (ic(0 : fc))] ^. Similarly, the prediction FIM J^'^ (a;(A; + l|A;)) is given by the 
inverse of the {n^ x Ux) right-lower block of [j^'^ (ic(0 : A;+l|fc))] ^. 

The problem we wish to solve is to compute the global information sub-matrix, denoted by J (x{k+l)) , 
as a function of the local FIMs j(') [x{k+l)) and local prediction FIMs jC) [x{k+l\k)) , for (1 < / < N). 
Note that J^^^ (^x{k)) can be updated sequentially using Eqs. (|421)-(|45]) where J{z^^\k + 1)) replaces 
J{z{k + 1)) in Eq. (l45l ). Proposition |2] derives a recursive formula for computing J^'^ (a;(A;+l|/c)) , i.e., 
the FIM for the local prediction distribution. 

Proposition 2. The sequence { J^'^ (a;(A;+ 1|A;)) } of the local prediction information sub-matrices for 
predicting state vectors x{k) at node I, for fl < / < N), follows the recursion 

j(')(a;(A:+l|A;)) = B^^{k) - D^\k){j^^\x{k)) + D^\k)y^ D^^ (k) (46) 

where j(')(a;(A;)) is given by Eq. dH]), D'^^{k), D^'^{k), and D'^^{k) are given by Eqs. and 

B^\k) = E{ - A^[^;+;j \ogP{x{k + 1)|^(A:))}, (47) 

The proof of Proposition |2] is included in Appendix |Dl Theorem |3] is our main result. It provides the 
exact recursive formula for computing the distributed FIM corresponding to the global estimation from 
the local FIMs j'^^\x{k)) and local prediction FIMs j'^^\x{k+l)). 

Theorem 3. The sequence {j(a;(A;))} of information sub-matrices corresponding to global estimates 
follows the recursion 

J{x{k+l)) = C^\k) - D^\k){j{x{k)) + D^\k)y^D^\k) (48) 
where D^^{k), D'^^{k), and D^'^{k) are given by Eqs. (|?i])-(|?5l) and 

N N 

C22(A;) = j;j«(^(A;+l))- j;j«(^(A;+l|A;))+E{-A^[^+;jlogP(^(A; + l)|:r(A;))}, (49) 
1=1 1=1 

where j('-*(a;(A;+l)) and J^''\x{k+l\k)) are defined in Prepositions \I}and^ respectively. 

The proof of Theorem |3] is included in Appendix O In |4 |, an approximate updating equation based on 
the information filter (an alternative form of the Kalman filter) is proposed for computing J{x{k + 1)) 
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at node I which is represented in our notation as follows 




Term J^'-\x{k+l)) is given by 

J^^\x{k + 1)) = D^\k) - D^\k)(^j'^^\x{k)) + D^\k)y^ D^\k). (51) 

Subtracting Eq. (46) from the above equation and rearranging, we get 

D^'^ik) = (^J^^\x{k+l))-J^^\xik+l\k))j +B^\k). (52) 

Substituting B'^'^{k) from Eq. (47) in Eq. (l52l ) and then substituting the resulting expression of D'^'^{k) 
in Eq. (|5TI) . we get 

J«(^r(A:+l))=(^jW(a;(fc+l))-J«(^(A:+l|A;))^+E{-A^j^+;j log P(^r(A; + 1)|^(A;)) } 

-D^\k){j'^^\x{k)) +D^\k)y^D^\k). (53) 

Substituting Eq (1531 ) in Eq. (l50l) results in the following approximated fused FIM at node I 

J{x{k + 1)) = C^\k) - D^\k){j'^^\x{k)) + D^\k)y^D^\k). (54) 

Note the differences between Eqs. (l54l ) and (l48l ). The second term in the right hand side of Eq. (l54l ) is based 
on the previous local FIM J^^\x{k)) at node / thus making it node-dependent, while the corresponding 
term in Eq. ( |48] ) is based on the overall FIM from the previous iteration. When the PCRLB is computed 
in a distributed manner, Eq. (l54l ) differs from one node to another and is, therefore, not conducive for 
deriving the overall PCRLB through consensus. To make Eq. dSOl ) node independent, our simulations also 
compare the proposed exact PCRLB with another approximate expression, which only includes the first 
two terms of C^^(A;) in Eq. (|49l) i.e., 

J{x{k+l)) ^ C^^k) ^ (j^^Hx{k + l))-J^^\x{k+l\k))\ . (55) 

1=1 ^ ^ 
The expectation terms in Eqs. (|43])-(|45]). ( |471 ). and (|49l ) can be further simplified for additive Gaussian 

noise, i.e., when ^(•) and C^'''(") normally distributed with zero mean and covariance matrices Q{k) 



and R^''\k), respectively. In this case, we have 

D'\k) = E{[V^^k)f^{k)]Q-\k)[V^^k)f''{k)f} (56) 

D'\k) = -E{[V^^,)f''{k)]Q-\k) = [D^\k)f (57) 

£>22(^)^Q-l(^)+E|[V^^^^,)/r(A; + l)]i?(0-^(A;+l)[V,(,+,)5W"(fe + l)]'^} (58) 

B'\k) = Q-\k), (59) 
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N N 

and C^\k) = jy^^\x{k+l))-Y/^^Hx{k+l\k)) + Q-\k). (60) 

1=1 1=1 

Note that Theorem [3] (Eqs. (l48])-(|49l)) provides a recursive framework for computing the distributed 
PCRLB. Further, the proposed distributed PCRLB can be implemented in a distributed fashion because 
Eq. ( |49l ) has only two summation terms involving local parameters. These terms can be computed in 
a distributed manner using the average consensus algorithms |[26l . Other terms in Eqs. (l48])-(|49l) are 
dependent only on the process model and can be derived locally. In cases (non-linear/non-Gaussian 
dynamic systems) where direct computation of D^^{k), D^'^{k), D'^^{k), D'^'^{k), B'^'^{k), and C'^'^{k) 
involves high-dimensional integrations, particle filters can alternatively be used to compute these terms. 

VI. Simulations: Bearing Only Target Tracking 



A distributed bearing-only tracking (BOT) application f391 is simulated to test the proposed CF/DPF. 
The BOT problem arises in a variety of nonlinear signal processing applications including radar surveil- 
lance, underwater submarine tracking in sonar, and robotics [39|. The objective is to design a practical 
filter capable of estimating the state kinematics x{k) = [X{k),Y{k),X{k),Y{k)\ (position [X, y] and 
velocity [X,y]) of the target from the bearing angle measurements and prior knowledge of the target's 
motion. The BOT is inherently a nonlinear application with nonlinearity incorporated either in the state 
dynamics or in the measurement model depending on the choice of the coordinate system used to 
formulate the problem. In this paper, we consider non-linear target kinematics with a non-Gaussian 
observation model. A clockwise coordinated turn kinematic motion model |[39l given by 



, p, sin(n(fc)AT) l-cosffi(fc)AT) 
g ^ l-cos(Q(fc)Ar) siii(n(fc)AT) 



with VL{k) = ^'"^ (61) 

^{x{k)Y^{Y(k)Y 



n{k) n(k) 
cos(0(A;)AT) -sm{n{k)AT) 

sm{n{k)AT) cos{n{k)AT) 

is considered with the manoeuvre acceleration parameter Am set to 1.08 x 10~^units/s^ ||39l . A sensor 
network of = 20 nodes with random geometric graph model in a square region of dimension (16 x 16) 
units is considered. Each sensor communicates only with its neighboring nodes within a connectivity 



radius of y/2 log{N)/N units. In addition, the network is assumed to be connected with each node 
linked to at least one other node in the network. Measurements are the target's bearings with respect to 
the platform of each node (referenced clockwise positive to the y-axis), i.e.. 
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where 

are the coordinates of node /. The observations are assumed to be corrupted by the 
non-Gaussian target glint noise fl31 modeled as a mixture model of two zero-mean Gaussians |45 |, one 
with a high probability of occurrence and small variance and the other with relatively a small probability 
of occurrence and high variance. The likelihood model at node /, for (1 < Z < N), is described as 

\x{k)) = (1 - e) X M{x- 0, (jj,, (k)) + e x M{x; 0, lOV^f,, (k)), (63) 

where e = 0.09 in the simulations. Furthermore, the observation noise is assumed to be state dependent 
such that the bearing noise variance (T^(,)(fc) at node I depends on the distance A^\k) between the 
observer and target. Based on P6l . the variance of the observation noise at node / is, therefore, given by 

0-2,,) {k) = O.OSrC)' {k) + 0.1150r(') {k) + 0.7405. (64) 

Due to state-dependent noise variance, we note that the signal to noise ratio (SNR) is time-varying and 
differs (within a range of — lOdB to 20dB) from one sensor node to the other depending on the location 
of the target. Averaged across all nodes and time, the mean SNR is 5.5dB. In our simulations, we 
chose to incorporate observations made at all nodes in the estimation, however, sensor selection based 
on the proposed distributed PCRLB can be used, instead, which will be considered as future work. Both 
centralized and distributed filters are initialized based on the procedure described in |[39ll . 
Simulation Results: The target starts from coordinates (3, 6) units The position of target the target 
{[X,Y]) in first three iterations are (2.6904,5.6209), (2.3932,5.2321), and (2.1098,4.8318). The initial 
course is set at —140° with the standard deviation of the process noise = 1.6 x 10^'^ unit. The number 
Ns of vector particles for centralized implementation is Ng = 10, 000. The number N^p and A'ff of vector 
particles used in each local filter and fusion filter is 500. The number of particles for the CF/DPF are 
selected to keep its computational complexity the same as that of the centralized implementation. To 
quantify the tracking performance of the proposed methods two scenarios are considered as follows. 
Scenario 1: accomplishes three goals. First, we compare the performance of the proposed CF/DPF versus 
the centralized implementation. The fusion filters used in the CF/DPF are allowed to converge between 
two consecutive iterations of the localized particle filters (i.e., we follow the timing subplot (a) of Fig. |2]). 
Second, we compare the impact of the three proposal distributions listed in Section III-E on the CF/DPF. 
The performance of the CF/DPF is computed for each of these proposal distributions using Monte Carlo 
simulations. Third, we plot the proposed PCRLB computed from the distributed configuration and compare 
it with its counterpart obtained from the centralized architecture that includes a fusion center. 

Fig. [3la) plots one realization of the target track and the estimated tracks obtained from: (i) The CF/DPF; 
(ii) the centralized implementation, and; (iii) a single node estimation (stand alone case ). In the CF/DPF, 
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Distributed vs Centralized PF Empirical CDF at iteration k = 5 Empirical CDF at iteration k = 22 
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-8 -6 -4 -2 2 4 6 8 

X - dimension 



(a) (b) 

Fig. 3. Scenario 1: (a) Actual target's track obtained from the centralized, CF/DPF and stand-alone algorithms. Here, the 
consensus algorithm is allowed to converge, (b) CDFs for the X-coordinate of the target estimated using the centralized and 
CF/DPF for fc = 5, 22. 



we used the Gaussian approximation of the optimal proposal distribution as the proposal distribution 
(Case 3 in Section III-E). The two estimates from the CF/DPF and the centralized implementation are 
fairly close to the true trajectory of the target so much so as that they overlap. The stand alone scenario 
based on running a particle filter at a single node (shown as the red circle in Fig. Oa)) fails to track 
the target. Fig. [3tb) plots the cumulative distribution function (CDF) for the X-coordinate of the target 
estimated using the centralized and CF/DPF implementations for iterations k = 5 and 22. We note that 
the two CDFs are close to each other. Fig. [3] illustrates the near-optimal nature of the CF/DPF. 

Fig. m compares the root mean square (RMS) error curves for the target's position. Based on a Monte- 
Carlo simulation of 100 runs, Fig. |4la) plots the RMS error curves for the estimated target's position 
via three CF/DPF implementations obtained using different proposals distributions. We observe that the 
SIR fusion filter performs the worst in this highly non-linear environment with non-Gaussian observation 
noise, while the outputs of the centralized and the other two distributed implementations are fairly close 
to each other and approach the PCRLB. Since the product fusion filter requires less computations, the 
simulations in Scenario 2 are based on the CF/DPF implementation using the product fusion filter. 
Second, in Fig. lUb), we compare the PCRLB obtained from the distributed and centralized architectures. 
The Jacobian terms Vx{k)f'^if^) Vx{k+i)9^''^^ {^ + ^)^ needed for the PCRLB, are computed using 
the procedure outlined in [39]. It is observed that the bound obtained from the proposed distributed 
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(a) (b) 

Fig. 4. Scenario 1: (a) Comparison of the RMS errors resulting from the centrahzed versus distributed implementations, (b) 
Comparison of the proposed exact PCRLB, its approximated value based on f4| (Eq. i54l ) computed at node 12 and 17, and 
approximation of Eq. i55\ . The PCRLBs computed using the centralized and distributed (proposed) expressions overlap so that 
they are virtually indistinguishable. 



computation of the PCRLB is more accurate and closer to the PCRLB computed from the centrahzed 
expressions. As expected, the proposed decentrahzed PCRLB overlaps with its centralized counterpart. 
Scenario 2: The second scenario models the timing subplot (c) of Fig. |2l The convergence of the fusion 
filter takes up to two iterations of the localized filters. The original fusion filter (Algorithm 1) is unable 
to converge within two consecutive iterations of the localized particle filters. Therefore, the lag between 
the fusion filters and the localized filters in the CF/DPF continues to increase exponentially. The modified 
fusion filter described in Algorithm |2] is implemented to limit the lag to two localized filter iterations. The 
target's track are shown in Fig.jSfa) for the centralized implementation, original and modified fusion filter. 
Fig. [5lb) shows the RMS error curves for the target's position including the RMS error resulting from 
Algorithm 1. Since consensus is not reached in Algorithm 1, therefore, the fusion estimate is different 
from one node to another. For Algorithm 1, Fig. |5tb) includes the RM error associated with a randomly 
selected node. In Fig.jSjb), Algorithm 1 performs poorly due to consensus not reached, while the modified 
fusion filter performs reasonably well. The performance of the modified fusion filter remains close to its 
centralized counterpart, therefore, it is capable of handling intermittent consensus steps. 
Scenario 3: In the third scenario, we consider a distributed mobile robot localization problem ifTSl . P8l 
based on angle-only measurements. This is a good benchmark since the underlying dynamics is non-linear 
with non-additive forcing terms resulting in a non-Gaussian transitional state model. The state vector of 
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(a) (b) 

Fig. 5. Scenario 2: (a) Actual target's track alongside estimated tracks obtained from the centralized and modified fusion filter. 
Here, the consensus algorithm converges after each 2 iterations of the local particle filters, (b) Comparison of the RMS errors 
resulting from the centralized, original fusion filter and modified fusion filter. 



the unicycle robot is defined by a; = [X, Y^O]^ , where (X, Y) is the 2D coordinate of the robot and 9 
is its orientation. The velocity and angular velocity are denoted by V{k) and W{k), respectively. The 
following discrete-time non-linear unicycle model HSl represents the state dynamics of the robot 

X{k+l) = X(A;) + ^^(sin(0(fc) + W^(A;)Ar) -sin(0(A;))), (65) 

Y{k + l) = Y{k) + ^!^^^{cos{e{k) + W{k)^T^ -cos{e{k))Y (66) 

and^(A;+l) = e{k) + W{k)/^T + ie^T, (67) 

where AT is the sampling time and is the orientation noise term. The design parameters are: AT = 1, 
a mean velocity of 30 cm/s with a standard deviation of 5 cm/s, and a mean angular velocity of 0.08 
rad/s with a standard deviation of 0.01 rad. Because of the presence of sine and cosine functions, the 
overall state dynamics in Eq. (I65l)-(l66l) are in effect perturbed by non-Gaussian terms. The observation 
model is similar to the one described for Scenario 1 with non-Gaussian and state-dependent observation 
noise. The robot starts at coordinates (3,5). Fig. |6] (a) shows one realization of the sensor placement 
along with the estimated robot's trajectories obtained from the proposed CF/DPF, centralized particle filter 
and distributed unscented Kalman filter (UKF) |fT5l implementations. We observe that both centraUzed 
particle filter and CF/DPF clearly follow the robot trajectory, while the distributed UKF deviates after 
a few initial iterations. Fig. [6] (b) plots the RMS errors obtained form a Monte-Carlo simulation of 100 
runs, which corroborate our earlier observation that the CF/DPF and the centralized particle filter provide 
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Fig. 6. (a) Robot trajectories estimated from the CF/DPF, centralized, and distributed UKF implementations, (b) RMS error 
plots for the three implementations. 



better estimates that are close to each other. The UKF produces a different result with the highest eiTor. 

VII. Conclusion 

In this paper, we propose a multi-rate framework referred to as the CF/DPF for distributed imple- 
mentation of the particle filter. In the proposed framework, two particle filters run at each sensor node. 
The first filter, referred to as the local filter, recursively runs the particle filter based only on the local 
observations. We introduce a second particle filter at each node, referred to as the fusion filter, which 
consistently assimilate local estimates into a global estimate by extracting new information. Our CF/DPF 
implementation allows the fusion filter to run at a rate different from that of the local filters. Achieving 
consensus between two successive iterations of the localized particle filter is no longer a requirement. The 
fusion filter and its consensus-step are now separated from the local filters, which enables the consensus 
step to converge without any time limitations. Another contribution of the paper is the derivation of the 
optimal posterior Cramer-Rao lower bound (PCRLB) for the distributed architecture based on a recursive 
procedure involving the local Fisher information matrices (FIM) of the distributed estimators. Numerical 
simulations verify the near-optimal performance of the CF/DPF. The CF/DPF estimates follows the 
centralized particle filter closely approaching the PCRLB at the signal to noise ratios that we tested. 

Appendix A 

Proof of Theorem |7] /|35|/.- Theorem [T] is obtained using: (i) The Markovian property of the state 
variables; (ii) Assuming that the local observations made at two sensor nodes conditioned on the state 
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variables are independent of each other, and; (iii) Using the Bayes' rule. Applying the Bayes' rule to 
Eq. ([T3]) . the posterior distribution can be represented as follows 

P{x{0:k)\z{l:k)) (xP{z{k)\x{k))P{x{0:k)\z{l:k-l)) . (68) 

Now, using the Markovian property of the state variables, Eq. (l68l) becomes 

P{x{0:k)\z{l:k)) (xP{z{k)\x{k)) x P{x{k)\x{k-l))P{x{0:k-l)\z{l:k-l)) . (69) 

Assuming that the local observations made at two sensor nodes conditioned on the state variables are 
independent of each other Eq. ( [69l ) becomes 

P{x{0:k)\z{l:k)) oc (^^P{z^^\k)\x{k))^ x P{x{k)\x{k-1)) P{x{0:k-l)\z{l:k-l)) . (70) 
Using the Bays' rule, the local likelihood function P (^z^^\k)\x{k)) at node I, for {1 < I < N) is 



Finally, the result (Eq. ([T3] )) is provided by substituting Eq. (TtTI ) in Eq. (iTOb . 



(71) 



Appendix B 

Proof of Theorem |2} Following the approach in the proof of Theorem [T] (Appendix |A]|, we first 
write the posterior density at iteration A;+m as 

ll^^^P {x{k+m)\z(^Hl:k+m)) 
nlli P {x{k+m)\z(^) {l:k+m-l)) 



P{x{0:k + m)\z{l:k+m)) oc J'"^^ , 77^^77777^- -f-:^P {x{0:k + m)\z{l:k + m-l)).(72) 



Then the last term is factorized as follows 

P {x{0:k + m)\z{l:k + m-l)) = P {x{k + m)\x{k + m-l)) P {x{0:k + m-l)\z{l:k + m-l)) . (73) 

As in Eq. (r72l ). we continue to expand P(a;(0:A;+m— l)|2;(l:A;+m — 1)) (i.e., the posterior distribution 
at iteration k+m — 1) all the way back to iteration + l to prove Eq. (l34l ). ■ 

Appendix C 

Proof of Theorem\3} Based on Eq. ([T3] ). term log(P(a;(0: A;+l)|z(l:/i;+l))) is given by 

N N 
\ogP{x{0:k+l)\z{l:k+l)) = ^log (^P{x{k+l)\z^^^l:k+l))^ - Z^log (^P{x{k + l)\z^^^ {l:k)) 

1=1 1=1 
+ log (P ix{k + l)\x{k)) ) + log {PixiO:k)\z{l:k))). 
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Expanding Eq. (|4T]) for J(a;(0: A;+l)), we get 



J{x{0:k+1)) = E< - 1 



x{k) 

x(0:k-l 

xlk+1) 



a; 

A^ 



x{k) 


.xik+1) 


x{0:k~l) 


^x{0:k-l) 


xlk) 


.x{k+l) 
^x{k) 


x{k) 


x{k) 


.x{k+l) 
^x{k+l) _ 


x{k+l) 



log P{x{0:k+l)\z{l:k+l)) \ (74) 



Substituting Eq. (1741) in Eq. (1741 ). it can be shown that 
J{x{0:k+1))-- 







(75) 



A^\k) A22(A:) + £)ii(A;) D^^{k) 
D^\k) C^^{k) 

where A^^{k), A^'^{k), A^^ik) and A^'^{k) are the same as for Eq. (SB; D^'^ik), D^'^{k) and D^^k) 
are the same as in Eqs. (I43l)-(l45l). and; C'^'^{k) is defined in Eq. ( |49l ). Block stands for a block of all 
zeros with the appropriate dimension. The information sub-matrix J {x{k+l)) can be calculated as the 
inverse of the right lower {n^ x rix) sub-matrix of [j(a;(0:A;+l))] ^ as follows 

A^^{k) A^'^{k) 



J{x{k+l)) = C22(A;)-ro D'^Hk)] 

A^Hk) A^\k) + C^\k) 
= C^^{k) - D^\k){j{x{k)) + D^\k)y^D^\k), 

obtained from the definition of J{x{k)) based on Eq. (|4T1) . 



-1 





D^'^{k) 



(76) 



Appendix D 

Proof of Proposition |2} The proof of Proposition |2] is similar to that of Theorem [3] using the 
Markovian property of the state variables and based on the factorization of the joint prediction distribution 
P{x{id:k+l)\z{l:k)) = P{x{k+l)\x{k))P{x{{):k)\z{l:k)) . ■ 
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